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combinations  thereof.  For  example,  in  cases  where  significant  diffraction  and  interference 
effects  require  •exact"  solutions,  finite  difference  techniques  have  received  widespread 
use.  v^Yet,  the  finite  difference  approach  is  well  known  for  its  cumbersome  computational 
denrfands  in  two  dimensions  and  almost  insurmountable  computational  demands  in  three  dimensions 
even  on  the  fastest  computers  available  today. 

The  heavy  computation  requirements  of  the  finite  difference  type  methods  are  created  by  the 
necessity  of  refining  the  numerical  grid  proportionately  to  the  wavelengths  of  interest  in 
all  spatial  directions,  including  regions  of  constant  material  properties.  For  problems 
involving  wave  propagation  in  irregular  layers  of  constant  material  properties  within  each 
layer,  however,  the  Boundary  Integral  Equation  (BIE)  approach  provides  a  more  concise  and 
efficient  formulation. 
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1.0  INTRODUCTION 


1.1  Review  of  Literature 


The  efficient  numerical  propagation  of  waves  in  complex 
three-dimensionally  varying  environments  has  been  a  problem  of 
considerable  geophysical  interest  over  the  past  few  years,  yet  has 
proven  extremely  difficult  even  for  the  case  of  acoustic  wave 
propagation  in  two  dimensional  structures.  This  is  primarily  due  to  the 
fact  that  although  the  differential  equations  of  motion  are  linear  in  the 
field  quantities  of  interest,  they  are  non-linear  in  terms  of  the 
boundary  conditions  for  most  realistic  structures.  This  fundamental 
nonlinearity  precludes  construction  of  the  solution  for  complex 
structures  by  superposition  of  the  solutions  for  simple  structures,  and 
forces  one  into  computationally  costly  schemes. 

Techniques  for  dealing  with  this  fundamental  nonlinearity  have  spanned 
the  range  from  the  crudest  classical  ray  tracing  approach  to  the 
computational-bound  finite  difference  type  methods.  However,  no  single 
technique  has  ever  proven  entirely  satisfactory  for  reasons  of  accuracy, 
completeness  of  solution,  generality  of  application,  cost  or  combinations 
thereof.  For  example,  in  cases  where  significant  diffraction  and 
interierence  effects  require  "exact"  solutions,  finite  difference 
techniques  have  received  widespread  use.  Yet,  the  finite  difference 
approach  is  well  known  for  its  cumbersome  computational  demands  in  two 
dimensions  and  almost  insurmountable  computational  demands  in  three 
dimensions  even  on  the  fastest  computers  available  today. 

The  heavy  computation  requirements  of  the  finite  difference  type 
methods  are  created  by  the  necessity  of  refining  the  numerical  grid 
proportionately  to  the  wavelengths  of  interest  in  all  spatial  directions, 
including  regions  of  constant  material  properties.  For  problems 
involving  wave  propagation  in  irregular  layers  of  constant  material 
properties  within  each  layer,  however,  the  Boundary  Integral  Equation 
(BIE)  approach  provides  a  more  concise  and  efficient  formulation. 


Basically,  the  BIE  formulation  takes  advantage  of  the  fact  that  the 
propagation  of  waves  through  a  region  of  constant  material  properties 
can  be  treated  analytically,  leaving  only  the  interactions  at  the 
bounding  surfaces  to  be  treated  numerically.  Rather  than  imposing  a 
grid  over  the  entire  volume,  the  BIE  method  only  requires  gridding  of 
the  interfaces  between  regions  of  constant  material  properties.  Not 
only  are  there  potential  savings  in  computational  effort  to  solve  a 
smaller  system  of  equations,  but  the  formulation  represents  a  concise 
treatment  of  the  pertinent  physics  involved.  By  virtue  of  this 
contraction  of  information,  the  smaller  matrices  in  the  BIE  approach  are 
much  denser  than  the  corresponding  matrices  in  the  finite  difference 
approach.  These  dense  matrices  are  typically  poorly  conditioned  and 
must  be  given  careful  consideration  during  implementation  of  matrix 
solution  techniques  to  avoid  numerical  instabilities. 

Various  techniques  have  appeared  in  the  literature  for  dealing  with  the 
dense  matrices  in  the  BlE  approach.  One  technique  involves 
introduction  of  a  Kirchhoff  approximation  into  the  BIE  formalism  (eg., 
Berryhill,  1979;  Scott  and  Helmberger,  1982;  Mellman,  et  al.,  1982).  In 
the  Kirchhoff  approximation  the  interaction  between  neighboring  points 
on  a  boundary  is  ignored  by  locally  approximating  the  boundary  at  each 
sample  point  by  the  tangent  plane  at  that  point  and  then  using  plane 
wave  reflection  and  transmission  coefficients  to  compute  the  unknown 
boundary  values.  Even  with  the  Kirchhoff  approximation,  one  is  still 
confronted  with  the  denseness  of  the  matrices  used  to  propagate  the 
boundary  values  forward  to  the  desired  positions.  Furthermore,  and  of 
utmost  importance,  is  the  fact  that  this  decoupling  of  neighboring 
boundary  points  in  the  Kirchhoff  approximation  precludes  simulation  of 
head  waves,  surface  waves,  most  diffraction  effects  and  any  other 
dynamic  effects  related  to  multiple  interactions  of  the  wavefield  with  a 
single  interface. 

A  time  domain  treatment  of  the  full  system  of  equations  has  been 
addressed  by  Cole  (1980)  hr  two-dimensional  acoustical  problems  in 
geophysics.  Cole's  approach  becomes  expensive  at  high  frequencies  or 
for  late  arriving  signals  as  the  product  of  the  frequency  step  times  the 
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time  step  must  be  less  than  about  10  to  maintain  stability.  More 
importantly,  the  formulation  does  not  handle  the  dense  matrices 
efficiently,  precluding  generalization  to  three-dimensional  elastic 
multilayered  problems.  Also,  it  is  difficult  to  include  realistic  material 
attenuation  and  to  suppress  late  arriving  spurious  reflections  off  the 
artificial  extremities  of  the  grid  using  a  time  domain  formulation. 

Ferguson  (1982)  studied  two-dimensional  elastic  problems  using  a 
frequency-domain  BIE  treatment  in  which  the  unknown  boundary  values 
are  expanded  in  a  series  of  plane  waves  with  unknown  amplitudes 
determined  by  performing  enormous  matrix  inversions  at  each 
frequency/wavenumber  pair  to  satisfy  the  boundary  conditions. 
Although  realistic  attenuation  is  included,  the  computational  cost  of 
Ferguson's  approach  is  at  least  an  order  of  magnitude  larger  than  finite 
difference  type  methods  and  provides  incorrect  results  for  problems 
involving  interfaces  with  slopes  exceeding  about  60  degrees. 

Schuster  (1984)  has  formulated  a  frequency  domain  BIE  approach  based 
on  first  solving  a  set  of  smaller  uncoupled  singular  boundary  integral 
equations  for  the  individual  primary  responses  of  each  interface  and 
then  coupling  them  together  by  successive  iterations  using  a  Neumann 
series  perturbation  treatment.  Schuster's  method  is  stable,  accurate, 
nicely  convergent  and  increases  in  cost  linearly  with  the  number  of 
layers,  yet  the  algorithm  still  requires  large  matrix  inversions  for  the 
individual  self-interaction  operators  preventing  a  cost-competitive 
alternative  to  finite  difference  methods. 

Apsel,  et  al.,  (1983)  formulated  a  frequency  domain  BIE  approach  in 
which  there  are  no  matrix  inversions,  realistic  attenuation  is  included, 
the  method  is  stable  and  accurate  and  the  algorithm  is  significantly 
more  cost-effective  than  finite  difference  type  algorithms.  A 
fundamental  ingredient  in  the  formulation  is  the  realization  that  the 
integrable  singularities  in  the  self-interaction  operators  along  each 
interface  have  the  same  convolutional  form  as  those  for  a  flat  reference 
plane.  Then  instead  of  perturbing  the  entire  primary  response  of  the 
individual  interfaces  as  in  Schuster's  approach  (1984),  the  exact  solu- 
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tion  including  all  kinematic  and  dynamic  effects  is  obtained 
iteratively  from  the  singular  self-interaction  responses  using  a  specially 
designed  perturbation  treatment  guaranteed  to  be  uniformly  and 
optimally  convergent.  All  matrix  inverse  operations  are  reduced  to 
simple  deconvolutional  operations,  which  are  efficiently  handled  using 
Fast  Fourier  Transform  algorithms. 


.•* 

•  •  *  •  ^ 
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1.2  WORK  COMPLETED  TO  DATE 

The  BIE  formulation  presented  in  Section  4.1  of  Apsel  et  al.  (1983)  has 
been  implemented  for  2-D  multilayered  acoustic  geologic  structures. 
There  were  two  significant  problems  with  the  original  implementation. 

First,  the  Neumann  series  iterative  procedure  exhibited  poor  and  often 
non-existent  convergence  for  models  departing  even  moderately  from  the 
flat  reference  planes.  The  second  problem  was  the  presence  of 
spurious  edge  reflections  for  models  with  interfaces  that  failed  to  return 
to  the  depths  of  the  reference  planes  near  the  horizontal  model 
extremes . 

To  address  the  convergence  problems,  it  was  necessary  to  make  four 
improvements  to  the  simple  Neumann  series  iterative  procedure.  The 
first  improvement  was  to  express  the  boundary  values  at  the  n-th 
iteration,  X^,  as  a  series  of  basis  vectors,  4>j,  with  unknown 
coefficients,  a.. 


1  a.  4'. 
i=1  '  ' 


-1  i-1  -1 

in  which  <l>j  =  ([C][A])  (C]{F}.  The  unknown  coefficients  at 

the  n-th  iteration  (a^ ,  0f2'  •••'  "p)  determined  by  minimizing  the 
residual  in  the  boundary  integral  equations  in  a  least-squared  error 
norm.  If  all  the  coefficients  were  determined  to  be  unity,  then  the 
expansion  in  Eq.  (1)  would  correspond  to  the  Neumann  series  solution. 
With  variable  coefficients  at  each  iteration,  the  method  is  guaranteed  to 
be  uniformly  convergent  in  the  absence  of  numerical  roundoff. 


The  second  improvement  was  to  orthogonalize  and  normalize  the  basis 
vectors  in  Eq.  (1)  in  order  to  more  rapidly  span  the  solution  space. 
This  resulted  in  approximately  a  20  percent  improvement  in  the  overall 
convergence  rate. 
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The  third  improvement  was  to  implement  an  automatic  restart  on  the 
iteration  loop  to:  eliminate  tne  potential  roundoff  problems;  reduce  the 
cost  of  the  least-squares  operation  at  ea<  h  iteration  by  limiting  the 

number  of  iterations  to  ten  per  pass;  and  increase  the  rate  of 
convergence  by  iterating  on  differences  of  boundary  values  from 
previous  passes  rather  than  directly  on  the  boundary  values.  This 

resulted  in  approximately  a  10  percent  improvement  in  the  overall 

convergence  rate. 

The  fourth  improvement  was  to  use  the  boundary  values  from  previous 

."I 

frequencies  to  achieve  a  better  starting  solution  than  [C]  {F}  at  the 

current  frequency.  This  resulted  in  approximately  a  20  percent 

reduction  in  the  number  of  iterations.  A  further  improvement  would  be 
possible  using  more  sophistic  ited  extrapolation  and  phase  unwrapping 
techniques  to  more  closely  predict  the  boundary  values  at  the  current 
frequencies  from  the  boundary  values  from  previous  frequencies. 

Even  though  these  four  improvements  provided  a  much  more  reliable 
algorithm,  the  convergence  was  still  far  too  slow  for  models  with 
moderate  or  large  perturbations  in  interface  depths  from  the  flat 
reference  planes. 

The  second  problem  area  related  to  spurious  edge  effects  was  addressed 
by  padding  the  models  by  at  least  ten  percent  at  both  horizontal 
extremes  and  applying  tapers  in  the  spatial  domain  to  suppress  edge 
reflections.  Wrap-around  events  in  the  spatial  domain  caused  by  the 
discrete  inverse  Fourier  transforms  in  the  wavenumber  deconvolution 
step  at  each  iteration  were  completely  suppressed  by  simple  padding  in 
the  wavenumber  domain.  Also,  potential  ringing  from  the  finite  Fourier 
transforms  was  suppressed  by  applying  tapers  at  large  wavenumbers  for 
the  inverse  transforms  and  in  the  spatial  taper  zones  for  the  forward 
transforms . 

This  combination  of  tapering  and  padding  was  very  effective  at 
eliminating  spurious  edge  effects  except  for  models  with  non-zero  relief 
from  the  reference  planes  near  the  edges  of  the  interfaces.  These 
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remaining  spurious  effects  and  the  inadequate  rate  of  convergence  are 
being  addressed  in  the  work  currently  in  progress  as  discussed  in 
Section  1.3.  The  original  approach  with  the  flat  reference  planes  is 
described  in  more  detail  in  Section  2.1  and  the  current  work  is 
described  in  more  detail  in  Section  2.2. 

Throughout  the  project,  rigorous  internal  and  external  validations  of 
the  algorithm  have  been  performed.  The  results  from  some  of  the  most 
important  validations  using  the  original  approach  are  presented  in 
Section  3.1.  Also,  some  preliminary  results  on  simple  models  for  AFGL 
are  presented  in  Section  3.2. 
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1.3  WORK  CURRENTLY  IN  PROGRESS 


The  outstanding  problems  discussed  in  Sec  tion  1.2  on  the  inadequate 
rate  of  convergence  and  the  spurious  eoge  effects  for  models  with 
interfaces  that  did  not  coincide  with  the  reference  planes  are  currently 
being  addressed.  Both  problems  were  directly  related  to  trying  to 
handle  general  models  with  large  perturbations  in  interface  depths  with 
respect  to  the  reference  planes  using  basiciilly  a  perturbation  approach. 
Using  the  flat  layer  deconvolutional  coefficients  to  precondition  the 
system  of  equations  did  not  improve  the  rate  of  convergence  for  models 
with  large  perturbations  and  did  not  solve  the  edge  effects  problem  for 
those  models  with  edge  perturbations  from  the  flat  reference  planes. 

To  address  both  problems,  the  new  formulation  has  eliminated  the 
dependence  on  the  reference  planes  to  suppress  edge  reflections  and  an 
improved  iterative  solution  technique  is  being  implemented  to  replace  the 
old  technique.  in  the  new  method,  the  models  are  padded  at  the 
horizontal  extremes  with  a  thin  absorption  zone  at  least  twenty  samples 
wide  in  which  the  forcing  functions  and  integration  quadrature 
coefficients  are  exponentially  tapered  to  zero  at  the  edges  to  suppress 
the  spurious  edge  reflections.  The  tapers  are  applied  in  the  frequency 
domain  by  prescribing  Q  values  that  are  tapered  to  neariy  zero  in  the 
absorption  zone.  This  is  similar  to  the  work  of  Cerjan,  et  al.  (1985) 
except  that  the  exponential  tapers  operate  only  on  the  amplitudes  and 
do  not  affect  the  phase  information. 


fm 


The  improved  iterative  solution  technique  is  an  asymmetric  conjugate 
direction  method  with  an  iteration  restart  capability  similar  to  the 
original  method  and  a  more  sophisticated  extrapolation  of  the  boundary 
values  from  previous  frequencies  to  achieve  a  closer  starting  guess  at 
the  next  frequency.  The  method  is  proving  to  be  stable  even  for  the 
most  complex  models  with  convergence  rates  on  the  order  of  the  square 
root  of  the  number  of  samples  per  interface.  The  largest  improvement 
is  the  more  uniform  convergence  rate  provided  by  more  optionally 
picking  new  search  directions,  whereas  the  rate  of  convergence  would 
slow  down  considerably  when  approaching  the  true  solution  in  the 
previous  method. 
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1.4  WORK  PLANNED  FOR  THE  NEXT  12  MONTHS 


The  first  task  is  to  fully  test  the  new  iterative  solution  technique  and 
new  edge  reflection  suppression  technique  in  the  2-D  acoustic  code  as 
described  in  Section  1.3.  This  will  involve  repeating  the  internal  and 
external  validation  exercises  using  the  upgraded  algorithm.  Once 
successful,  various  2-D  acoustic  simulations  of  interest  to  AFGL  will  be 
performed . 


The  next  phase  of  the  project  would  then  be  to  extend  the  algorithm  to 
the  3-D  acoustic  case  and  perform  more  validation  exercises  and 
simulations  for  3-D  cases.  After  the  3-D  acoustic  algorithm  is  complete, 
development  will  begin  on  the  3-D  elastic  algorithm. 
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2.0  METHODOLOGY 


2.1  ORIGINAL  BINTEQ  FORMULATION 

The  boundary  integral  equations  describing  complete  wave  propagation 
through  arbitrary  three-dimensional  elastic  multilayered  media  are 
derived  in  two  steps.  First,  the  known  characterization  of  wave 
propagation  within  a  single  irregular  layer  is  written  in  terms  of 
integral  representations  involving  the  full  space  Green's  functions  with 
properties  of  that  layer.  Second,  the  interaction  of  the  wavefield  is 
simultaneously  imposed  at  all  boundaries  to  satisfy  all  boundary  and 
continuity  conditions  leading  to  a  system  of  Fredholm  integral  equations 
of  the  second  kind  for  the  unknown  boundary  values.  Once  this 
system  of  equations  is  solved  for  the  unknown  boundary  values,  the 
wavefield  may  be  propagated  from  the  boundaries  to  all  receiver 
positions  of  interest  within  a  given  layer  using  the  integral 
representations  of  the  first  -.tep.  The  formulation  for  the  2-D  and  3-D 
acoustic  cases  is  analogous  to  the  3-D  elastic  case  and  will  not  be 
repeated  here. 

The  model  geometry  for  the  wave  propagation  problem  solved  in  the  BIE 
formulation  is  depicted  in  Figure  1  by  N  irregular  layers  overlying  a 
semi-infinite  half-space.  The  layers  are  allowed  to  pinchout  but  not  to 
cross  in  this  formulation.  Each  layer  is  characterized  by  constant 
shear  and  compressional  wave  velocities  and  constant  densities. 
Realistic  material  attenuation  is  introduced  by  allowing  the  velocities  to 
be  complex.  Wave  propagation  within  a  given  layer  is  expressed  in 
terms  of  the  Green's  functions  for  a  full-space  with  the  properties  of 
that  layer.  The  formulation  is  not  restricted  to  constant  material 
properties  within  a  given  layer,  although  the  Green's  functions  for  that 
case  are  quite  simple.  The  formulation  is  carried  out  for  the  full 
elastic  case  and  the  corresponding  acoustic  formulation  is  obtainable 
from  the  derived  equations  by  replacing  the  vector  equations  with 
scalar  equations. 
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The  first  step  in  the  formulation  is  to  write  expressions  for  the 
displacement  field  within  a  single  layer  without  consideration  of  the 
boundary  interaction.  In  layer  SL  (£=1,2, . . .  ,N+1),  the  displacement 
vector  must  satisfy  the  homogeneous  (£i^s)  or  inhomogeneous  (£=s) 
equations  of  motion  (depending  on  whether  or  not  the  source  layer  s  is 
the  same  as  layer  £)  for  a  full-space  with  properties  of  layer  £.  The 
Representation  Theorem  of  Elastodynamics  (see,  for  example,  deHoop, 
1958)  provides  an  expression  for  the  displacement  vector  located  any¬ 
where  within  volume  containing  layer  £  in  terms  of  integrals  of  the 
displacement  and  traction  field  over  the  bounding  surface  of  volume  \J ^ 
times  the  corresponding  Green's  functions  for  a  full-space  with 
properties  of  that  layer.  The  i-component  of  displacement  at  location 
can  then  be  written  in  the  frequency  domain  using  the 
Representation  Theorem  for  a  volume  bounded  by  layer  interfaces 


£,->  ,  ,,£,■♦  s 

e  (x^  U,(xp 


as(y,) 


’^£+1 


''  j  '•^£+1' 


r,(x^,2^)f.(z^)j 


dV(z^) 


(2) 


(i,j=l ,2,3) 


in  which  the  summation  over  repeated  indices  is  understood, 
frequency  arguments  have  been  omitted  for  brevity  and 


the 


an  integration  point  on  bounding  surface  S^; 


the  j-component  of  the  full-space  Green's^ 
function  displacement  vector  at  location  y 

*  lU 

on  surface  S  due  to  poi^t  force  in  the 
i-direction  at  location  x.  with  properties  of 
layer  £; 


SGI-R-85-120 


H. .  (x.,y  ) 

ji 


U^Cy  ) 


the  j-component  of  the  corresponding  Green's 

function  traction  vector  formed  from  the 

inner-product  of  the  kj-componenj^  of  the  ^ 

Green's  function  stress  tensor  Gj^.^  at  location  y^ 

on  surface  S  with  the  k-componeni^of  the 
,  m  ,  m  •  ^ 

unit  upward  normal  v,  at  point  y  (summing 

over  k=l ,2,3) ;  ™ 

^he  j-component  of  displacement  at  location 
y  on  surface  S  ; 

»  m  m  " 


T'  (y  ) 


the  j -component  of  the  corresponding  traction 
at  location  y  on  surface  S  formed  from  the 
inner-product  of  the  kj -component  of  the 
stress  tensor  with  the  k-compon^nt  of  the 
unit  upward  normal  v.  at  point  y  (summing 
over  k=l,2,3); 


=  the  j-comjjonent  of  the  source  function  at 

location  anywhere  in  layer  SL  (assuming  the 
source  is  a  Delta-function  in  space,  then  the 
volume  integral  reduces_^to  the  evaluation  of 
the  integrand  at  point  z^) ; 


is 


0,  if  l^s 
1,  if  l=s 


,s  =  source  layer  number; 


e  (x^) 


1 ,  if  inside  layer  £ 

if  X-  on  surface  bounding  layer  £ 

_ 1 _  A 


0, 


if  Xj^outside  layer  £. 


In  Eq.  (2),  the  layer  comprising  volume  is  assumed  to  extend  to 
infinity  at  the  horizontal  extremes  to  eliminate  the  surface  integrals 
along  those  portions  of  the  surface  bounding  volume  and  the  nega¬ 
tive  sign  for  the  integral  over  surface  -  is  associated  with  using  the 
^£+1  •♦£  ^ 

upward  normal  v  =  -v  in  the  definition  of  the  traction  components. 

Once  the  boundary  values  for  u’T\y  )  and  T^(y  )  are  determined  for 

J  m  J  fTi 

bounding  interfaces  and  EQ-  (T)  can  then  be  used  to  obtain 

the  displacement  field  at  any  point  x^  within  layer  £.  Expressions  for 
the  full-space  Green's  functions  with  constant  material  properties  are 
given  in  Appendix  A  of  Apsel,  et  al.  (1983)  for  two  and  three 
dimensional  wave  propagation  in  elastic  as  well  as  acoustic  media.  This 
completes  the  propagation  step  of  the  BIE  formulation  and  what  remains 
is  to  impose  the  boundary  interaction  coupling. 


The  boundary  interaction  coupling  requires  simultaneous  satisfaction  of 
a  tractionless  free  surface  (interface  1)  and  continuous  displacements 
and  tractions  across  each  layer  interface  (2,3, . . .  ,N+1).  The  coupled 
boundary  integral  equations  arising  from  the  zero  traction  conditions 
along  the  free  surface  are  obtained  by  evaluating  Eq.  (1)  in  volume 
(layer  1 )  at  a  discrete  number  (q,)  of  observation  points  x-  along 
surface  and  imposing  the  zero  traction  condition,  =  0/ 

(j=1,2,3),  for  all  quadrature  points  y^  on  surface  .  This  leads  to  a 
simultaneous  set  of  3q-  Fredholm  integral  equations  of  the  second  kind 
for  the  same  number  of  unknown  displacement  boundary  values  U., 
j=1,2,3,  on  surface  S^,  which  are  coupled  to  the  unknown  boundary 
values  on  surface  S2  through  the  integral  over  surface  Sg- 


T  (x^, 


The  coupled  boundary  integral  equations  arising  from  the  continuity 

conditions  across  each  layer  interface,  (£=2,3, . . .  ,N+1 ),  are  obtained 

by  evaluating  Eq.  (2)  in  volumes  V^_.j  and  V ^  (layers  £-1  and  £)  at  a 

discrete  number  (q£)  of  observations  x^  along  common  surface  and 

imposing  the  natural  boundary  conditions  of  continuous  displacements 

and  tractions,  l/  ^x^)  =  U^(x£)  and  T^ 

j=1,2,3,  for  all  quadrature  points  y^  on  surface  S^.  This  leads  to  a 

simultaneous  set  of  Gq^  Fredholm  integral  « quations  of  the  second  kind 

for  the  same  number  of  unknown  displacement  and  traction  boundary 

values,  1/  and  T^,  j=l,2,3,  on  surface  S,,  which  are  coupled  to  the 

unknown  boundary  values  on  surfaces  S£^.^  through  the 

integrals  over  surfaces  ,  and  S*  respectively.  Note  that  when 

^-1  -♦  ^ 

£=2,  the  integrals  involving  T:  (Vg.-j)  S'"®  identically  zero  because  of 
the  tractionless  free  surface  conditions.  Also,  note  that  when  £=N+1 , 
the  integrals  over  surface  vanish  by  virtue  of  the  radiation 

conditions  implicit  in 
semi-infinite  space. 


£+1 

the  Green's  functions  for  the  underlying 


When  the  entire  discrete  set  of  boundary  and  continuity  conditions  is 
simultaneously  imposed,  one  obtains  a  coupled  system  of  singular 
Fredholm  integral  equations  of  the  second  kind  for  the  unknown 
boundary  values  along  all  the  interfaces,  which  can  be  written  in  matrix 
notation  as: 
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[G]{T} 


[H]{U}  +  {F} 


(3) 


ll2](U}  = 

in  which  [l^]  is  a  di-diagonal  matrix  consisting  of  the  =  1/2  factors 
obtained  when  specializing  Eq.  (2)  to  points  on  the  interfaces^  [G]  and 
[H]  are  the  block  tri-diagonal  displacement  and  traction  Green's 
function  matrices,  respectively;  {F}  is  the  forcing  vector  consisting  of 
the  direct  source  contributions  at  nodes  only  on  the  interfaces  bounding 
the  source;  and  {U}  and  {T}  are  the  unknown  displacement  and  traction 
boundary  value  vectors,  respectively,  at  all  nodes  in  the  model.  The 
singularities  occur  in  [G]  and  [H]  when  quadrature  point  y^ 
approaches  observation  point  in  the  self-interaction  integrals  along 
each  interface  and  in  the  propagation  integrals  between  adjacent 
interfaces  for  the  special  case  of  a  layer  pinchout. 

The  original  BINTEQ  solution  to  Eq.  3  was  formulated  to  meet  four 
objectives: 

1)  optimize  computational  speed; 

2)  minimize  memory  requirements  for  efficient  execution  in 
array  processors  and/or  multi-user  environments; 

3)  suppress  spurious  edge  effects  from  the  horizontal 
finiteness  of  the  numerical  grid; 

4)  generate  accurate  complete  sclutions  including  all  possible 
kinematic  and  dynamic  effects; 


To  accommodate  these  objectives,  the  BINTEQ  formulation  proceeded  by 
recognizing  that  the  singular  self-inferaction  elements  of  matrices  [G] 
and  [H]  are  identical  in  the  limit  to  those  for  a  flat  reference  plane 
with  the  same  upper  and  lower  material  properties  as  the  interface 
being  considered.  The  rows  of  the  self-interaction  matrices  for  layer 
interfaces  are  simple  convolutional  operators  and  are  ideal  for 
preconditioning  the  original  matrix  equation.  Therefore,  by  subtracting 
the  self-interaction  matrices  for  flat  planes  referenced  to  each  interface 
from  both  sides  of  Eq.  (3),  the  system  of  equations  can  be  rewritten 


(4) 


[C){(t)}j^^  =  [A]{<^}.  +  {F} 

in  which  [C]  is  a  block  diagonal  matrix  containing  the  1/2  factors  from 
[I2]  and  the  singular  convolutional  coefficients  from  the  flat 
self-interaction  matrices;  [A]  is  the  combined  perturbation  Green's 
function  matrix  obtained  by  subtracting  the  flat  self-interaction  matrices 
from  [G]  and  [H];  and  {$}  is  the  combined  vector  of  unknown 
displacement  and  traction  boundary  values  {U}  and  {T}. 

For  models  with  small  perturbations  from  the  flat  reference  planes,  the 
preconditioned  system  of  equations  in  Eq.  (4)  is  more  numerically 
tractable  and  is  well  suited  to  satisfy  all  four  objectives  simultaneously. 
The  interaction  singularities  are  analytically  integrable  and  appear  only 
in  the  convolutional  matrix  [C]  and  the  pinchout  singularities  (if  any 
exist)  are  tIso  analytically  integrable  and  appear  only  in  perturbation 
matrix  [A],  If  it  could  be  guaranteed  that  the  norm  of  the  [A]  matrix 
is  always  less  than  the  norm  of  the  [C]  matrix,  then  Eq.  (4)  could  be 
solved  as  accurately  as  desired  using  the  following  iterative  Neumann 
series  solution  technique  with  =  0  to  initialize  the  series:  first  the 

solution  from  iteration  i,  is  recursively  passed  through  the 

right-hand-side  of  Eq.  (4);  then  the  right-hand-side  is  transformed 
into  the  horizontal  wavenumber  domain  using  a  discrete  FFT  algorithm; 
then  an  updated  solution  is  immediately  obtained  by  applying  the 

deconvolulional  coefficients  of  the  [C]  matrix  in  the  wavenumber 
domain;  and  then  the  updated  solution  is  transformed  back  into  the 
spatial  domain  for  the  next  it^mation. 

This  procedure  would  satisfy  the  first  objective  by:  (a)  eliminating  all 
matrix  inversion  operations;  (b)  saving  the  nonzero  perturbation 
submatrices  of  [A]  in  the  spatial  domain  and  the  nonzero 
deconvolutional  coefficients  of  |C]  in  the  wavenumber  domain  for 
recursive  iterations  and  multiple  sources;  (c)  reducing  the  number  of 
iterations  of  the  precondition  system  relative  to  the  original  system;  and 
(d)  using  an  array  processor  to  rapidly  calculate  the  nonzero  elements 
of  (A|  and  (Cj,  process  all  the  required  FFTs  and  perform  all  the 
required  complex  matrix/vector  multiplies. 


This  procedure  also  satisfies  the  second  objective  by:  (a)  saving 
nontrivial  perturbation  submatrices  and  deconvolutional  coefficients  on 
disk  if  memory  is  insufficient;  and  (b)  the  largest  in-core  memory 
requirement  is  governed  merely  by  a  single  complex  multiplication  of  a 
submatrix  times  a  bounoary  value  subvector. 

To  understand  how  the  third  objective  is  satisfied,  it  is  instructive  to 
consider  the  origins  of  the  three  possible  types  of  spurious  edge 
effects.  First,  edge  reflections  from  the  deconvolutional  operation  on 
the  forcing  vector  {Fj  during  the  first  iteration  are  possible  for  forcing 
vectors  without  compact  support,  which  would  usually  be  the  case 
except  possibly  for  structures  with  significant  amounts  of  material 
attenuation  in  the  source  layer  (i.e.,  low  Q  values).  Second,  edge 
reflections  are  similarly  possible  when  updating  the  right-hand-side  of 
Eq.  (4)  if  the  perturbed  Green's  function  integration  operators  in  [A] 
do  not  have  compact  support.  Third,  spatial  wraparound  effects  are 
possible  if  the  convolutional  coefficients  are  not  sufficiently  padded  with 
zeroes.  It  should  also  be  pointed  out  that  for  interfaces  which  return 
to  their  respective  flat  reference  planes  at  both  horizontal  extremes, 
the  second  type  of  edge  reflections  would  require  less  care  than  the 
first  type.  Therewith,  the  second  objective  is  correspondingly  satisfied 
by:  (a)  extending  the  model  somewhat  at  the  horizontal  extremes  with 

and  tapering  the  forcing  functions  and  perturbed  integration  operations 
to  zero;  (b)  tapering  the  convolutional  coefficients;  and  (c)  padding 
the  convolutional  coefficients  and  the  right-hand-sides  with  zeroes  out 
to  twice  the  model  size  to  totally  prevent  the  circular  deconvolution 
process  in  the  wavenumber  domain  from  wrapping  any  arrivals  back  into 
the  model. 

As  mentioned  previously,  if  one  could  guarantee  that  the  norm  of  [A] 
be  less  than  the  norm  of  [C],  then  the  iterative  Neumann  series 
solution  technique  would  converge  rapidly  and  uniformly  to  the  exact 
solution  and  the  fourth  objective  would  be  met  automatically.  However, 
this  cannot  always  be  guaranteed  especially  for  large  model 
perturbations  away  from  the  flat  reference  planes.  In  Schuster's 
iterative  BIE  solution  technique,  the  Neumann  series  is  guaranteed  to 
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be  uniformly  convergent  because  the  full  interaction  submatrices  are 
inverted,  leaving  only  l  ,e  coupling  between  interfaces  to  control  the 
rate  of  convergence.  As  will  be  seen  shortly,  however,  Schuster's 
approach  is  extremely  inefficient  because  of  having  to  perform  the 
matrix  inversions  and  furthermore,  there  are  alternatives  to  the 
Neumann  series  expansion  which  are  guaranteed  to  be  uniformly 
convergent.  The  alternative  adopted  for  the  original  BINTEQ  formulation 
is  to  expand  the  unknown  boundary  values  in  a  series  of  basis  vectors 
with  unknown  coefficients  as  shown  in  Eq.  (1)  in  Section  1.2.  Each 
basis  vector  is  generated  recursively  as  described  above  for  the 
Neumann  method  and  is  made  orthogonal  to  all  previous  basis  vectors 
using  a  modified  Gram-Schmidt  orthogonalization  procedure.  The 
unknown  basis  vector  coefficients  are  determined  at  a  given  iteration  to 
satisfy  the  boundary  and  continuity  conditions  implicit  in  Eq.  (4)  in  a 
least-squared  error  norm.  To  avoid  numerical  roundoff  problems,  an 
automatic  restart  on  the  iteration  loop  is  required  whenever  the 
condition  number  of  the  least-squares  system  for  the  unknown 
coefficients  exceeds  single  precision  accuracy.  The  iterated  solution 
contains  all  the  possible  arrivals  (e.g.,  direct  waves,  multiples, 
converted  phases,  head  waves,  diffractions  and  surface  waves).  Once 
all  the  boundary  values  are  determined  at  a  given  frequency,  the  field 
at  any  location  x^  within  any  layer  £  may  be  obtained  by  evaluating  Eq. 
(2).  Time  domain  response  would  be  obtained  through  discrete  Fourier 
synthesis . 

The  main  problem  with  this  procedure  is  the  requirement  that  the 
perturbations  of  the  irregular  interfaces  from  the  flat  reference  planes 
be  small  at  the  center  of  the  model  and  zero  at  the  edges.  Otherwise, 
the  preconditioning  from  subtracting  the  flat  self-interaction  matrices 
would  not  improve  the  rate  of  convergence  and  the  simple  tapering 
described  above  would  be  insufficient  to  suppress  all  the  spurious  edge 
effects.  Because  the  basis  vectors  are  recursively  dependent  on  the  fif* 

-I 

layer  self-interaction  solution  (fC)  {F}),  the  more  iterations  required 
to  satisfy  the  desired  error  tolerance,  the  slower  the  rate  of 
convergence.  Therefore  the  rate  of  convergence  is  model  dependent  and 
is  too  slow  for  practical  applications  using  this  original  procedure. 
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2.2  IMPROVED  BINTEQ  FORMULATION 

The  preconditioning  of  the  system  of  equations  in  Eq.  (4)  is  physically 
motivated  and  useful  only  for  small  perturbation  problems  with  respect 
to  the  flat  reference  planes.  To  solve  the  more  desirable  general 

layered  cases  more  efficiently,  it  was  necessary  to  choose  a  more 

reliable  preconditioning  technique  and  an  iterative  scheme  that  more 
rapidly  spans  the  solution  space. 

First,  the  flat  layer  reference  planes  are  being  eliminated  and  a  more 
robust  technique  is  being  implemented  to  suppress  the  spurious  edge 
reflections.  The  layers  are  padded  at  the  horizontal  extremes  with  thin 
absorption  zones  in  which  the  Q  values  for  the  layer  are  smoothly 
tapered  to  a  small  waterlevel  value  at  the  edge  of  the  model.  This 

introduces  an  extra  exponential  decay  into  the  forcing  functions  and 
Green's  function  integration  operators  which  gradually  reduces  the 
amplitudes  of  the  boundary  values  in  the  absorption  zone.  The 

waterlevel  value  is  set  at  each  frequency  such  that  any  spurious  edge 
reflections  are  Loo  small  to  contaminate  the  real  signal. 


The  implementation  of  the  absorption  zones  is  portrayed  in  Figure  2  for 
observation  point  in  the  absorption  zone  and  source  point  in  the 
original  unpadded  model.  Defining  to  be  the  width  of  the  absorption 
zone  (typically  the  greater  of  20  samples  or  10  percent  of  the  model 
width  in  that  direction),  the  Q-values  at  position  S'p.  in  the  absorption 
zone  are  given  by 


(l-w,) 


Y 


cos  n  Jlj 
2  A 


W, 


(5) 


in  which  "W  "  is  the  waterlevel  factor,  s  controls  the  power  of  the 
decay  and  Y^  signifies  the  horizontal  starting  position  of  one  of  the 
absorption  zones  or  X  ,  if  X  is  within  the  absorption  zone. 
Therefore,  the  Q-values  are  smoothly  tapered  from  a  value  of  at  the 
start  of  the  absorption  zone  to  WQ^  at  the  end  of  the  absorption  zone. 


Absorbtion 

Zbne 


Absorbtion 

Zone 


%  % 


inlerface^ 


Layer  i 

®/  •  Pi  >  Qi 

interface 


Layer, 


Aj  — 


Figure  2 

Geometry  and  notation  used  for  horizontal  absorbtion  zones. 
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The  differential  amplitude  reduction  factors,  RF,  for  the  Green's 
functions  are  then  given  by  multiplying  the  N  reduction  factors  for  the 
N  cells  between  and  Y^: 


for  which  lu  is  the  angular  frequency  and  Ar  is  the  slant  distance  along 
each  of  the  cells  contributing  to  the  reduction  factor.  The  waterlevel 
factor,  is  computed  to  make  the  reduction  less  than  or  equal  to  0.1 

for  the  last  cell's  contribution; 

(0.1)  ^  , 

wax 

If  AX  is  different  from  AY  for  3-0  models,  then  would  be  calculated 
separately  in  the  two  horizontal  directions.  The  power  s^^  is  computed 
such  that  Q^.  =  0.1  '^he  mid-point  of  the  absorption  zone  in 

order  to  provide  gradual  reduction  factors  throughout  the  absorption 
zone. 

This  procedure  to  suppress  edge  reflections  is  similar  to  the  work  of 
Cerjan,  et  al.  (1985),  with  two  exceptions.  First,  the  constants  in  the 
exponential  reduction  factors  are  based  on  physical  quantities  in  the 
present  approach  and  do  not  alter  the  phase  information.  Second, 
there  is  no  need  for  an  absorbing  zone  at  the  bottom  of  the  models 
because  the  radiation  conditions  are  handled  exactly  by  the  Green's 
functions  for  the  underlying  half-space  layer.  This  new  procedure  is 
in  the  process  of  being  tested  and  some  of  the  constants  may  need  to 
be  altered  for  optimal  suppression  of  the  spurious  edge  reflections. 
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Along  with  the  new  edge  reflection  suppression  procedure,  a  new 
iterative  solution  technique  is  being  implemented  and  tested.  The 
largest  problems  with  the  original  iterative  solution  technique  were; 
(1)  using  basically  a  perturbation  method  on  models  with  moderate  or 
large  perturbations  in  the  interface  depths  with  respect  to  the  flat 
reference  planes;  (2)  less  than  optimal  preconditioning  of  the  matrix 
equations  for  these  non-perturbation  problems;  and  (3)  no  facility  to 
pick  optional  search  directions  for  the  successive  iterations.  In  the 
new  method,  a  conjugate  gradient  iterative  solution  technique  is  used 
instead  and  the  singular  diagonal  terms  are  used  to  precondition  the 
system  of  equations. 

Numerically,  the  problem  with  the  original  method  was  that  the  search 
directions  became  too  similar  at  successive  iterations  causing 
prohibitively  slow  convergence  in  many  cases.  This  is  the  same  type  of 
problem  encountered  in  using  the  method  of  steepest  descent  where 
minimization  in  the  gradient  direction  causes  convergence  back  and 
forth  across  the  valley  rather  than  more  directly  down  the  valley.  The 
conjugate  gradient  (CG)  method  provides  a  framework  for  picking  the 
search  directions  to  minimize  the  residuals  more  rapidly  while  still 
guaranteeing  uniform  convergence. 

To  derive  the  CG  solution,  Eq.  (3)  is  recast  into  the  simple  form; 

A  X  =  B  (8) 

in  which  A  is  a  general,  asymmetric  complex  matrix  containing  the 
Green's  function  integ^'ation  submatrices  from  [G]  and  [H]  and  the 
£-factors  from  [I2],  B  is  the  forcing  vector  {Fj,  and  X  is  the  unknown 
boundary  value  vector.  The  standard  CG  method  is  for  symmetric 
positive  definite  matrix  equations.  Before  deriving  the  generalization  to 
the  asymmetric  case,  the  basic  formulae  and  properties  of  the  CG 
solution  will  be  discussed  for  the  normal  equations; 

A^AX  =  A^B  (9) 
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in  which  superscript  T  indicates  the  complex-conjugate  transpose 
matrix.  A  good  reference  for  the  symmetric  CG  solution  to  Eq.  (9)  can 
be  found  in  Chapter  10  of  Golub  and  van  Loan  (1983).  The  basic 
approach  is  to  span  the  solution  X  by  a  set  of  mutually  A-orthogonal 
search  directions  P^,  (n=1,2,...) 


X  =  X,+aP  ,  X  =  best  estimate, 
n  n-1  n  n  '  o 


(10) 


with  the  corresponding  residual  vector  r^ 


=  A^(B-AX^ 


)  given  by 


r 

n 


n-1 


a  A^AP 
n  n 


r  =  A^(B-AX  ) 
o  o' 


(11) 


To  directly  minimize  the  residual  vector  in  Eq.  (11),  the  coefficients 
are  found  by  requiring  that  (Pj^/  fp)  =  0  to  give 


“  -  ^n-1)  .  (''n-1,^n-1)  (12) 

^  ■  (AP  ,AP  )  (AP  ,AP  ) 

^  n'  n  n'  n 

in  which  the  second  form  is  derived  by  using  Eq.  (14)  and  induction 
arguments  and  the  notation  (X,Y)  denotes  the  inner  product  X  Y. 
What  remains  is  how  to  define  the  optimum  search  directions  that  satisfy 
the  A-orthogonality  condition: 


(AP.,  AP.)  =  0  for  i  ^  j  (13) 

To  reduce  the  residuals  as  rapidly  as  possible,  it  is  desirable  to  choose 

the  search  directions  P  to  be  the  closest  vectors  to  r  ,  that  still 

n  n-1 

satisfy  Eq.  (13).  With  no  loss  in  generality,  the  search  directions  can 
be  written  recursively  as: 


P 

n 


n-1 


+ 


P 


n 


P 

n- 


1 


(14) 
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Applying  Eq.  (73)  with  i  =  n,  introducing  from  Eq.  (14)  and  A  AP. 

from  Eq.  (11),  using  induction  arguments  for  j<n  and  the  orthogonality 

of  the  residual  vectors  (r.,n)  =  0  for  i  /  j  directly  gives  the 

coefficients  fi  : 

n 


(’'n-l,  ’’n-D 


(n  =  2,3,...). 


To  optimize  the  rate  of  convergence,  the  system  of  equations  in  Eq. 
(11)  is  preconditioned  by  the  diagonal  elements  and  the  initial  estimate 
is  derived  by  extrapolating  the  amplitude  and  phase  information  from 
the  solution  at  previous  frequencies. 


Note  that  this  CG  algorithm  for  the  normal  equations  requires  two 
matrix-vector  multiplications  for  each  iteration  because  of  the  A^AP 

n 

term.  An  additional  drawback  of  solving  the  normal  equations  is  that 
the  rate  of  convergence  is  governed  by  the  square  of  the  condition 
number  of  the  A-matri  <  instead  of  just  the  condition  number  with  an 
asymmetric  algorithm. 


The  derivation  of  the  asymmetric  CG  method  is  similar  to  the  symmetric 

case  in  Eqs.  (10)  through  (15)  with  a  few  basic  changes.  As  before, 

the  solution  is  expanded  in  a  set  of  A-orthogonal  search  directions  P^, 

(n=1,2,...)  with  coefficients  a  : 

'  '  n 


X  =  X  +  a  P 
n  n-1  n  n 


X^  =  best  estimate. 


with  residual  vector  r^  =  B-AX^  given  by 


=  r  ^  -  a  AP 
n-1  n  n 


r  =  B-AX 
o  o 
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In  this  case  the  coefficients  a  are  determined  to  minimize  the  residual 

n 

vector  in  Eq.  (17)  in  a  least-squared  error  sense  to  give: 


“n  “  (*'n,  ''n-D 

(APr,  AP„) 


Without  symmetry  arguments,  the  search  directions,  P^,  depend  on  all 
previous  search  directions  to  satisfy  Eq.  (13): 


P  =  r  +  I  p.  P.  ,  P,  =  r  .  (19) 

n  n-1  .  ^in  i  '  1  o 

i=1 

Now,  applying  Eq.  (13)  with  i=n  and  j  =  1,2,..., n-1  gives  n-1 
decoupled  equations  for  the  n-1  coefficients  p.^  at  iteration  n: 


o  ,Ar  .  AP.s 

Pin  =  < _ _ L) 

(AP.,  APj) 


i=1,2, . . . ,n-1 


This  asymmetric  procedure  has  the  distinct  advantage  of  only  one 
matrix-vector  multiplication  per  iteration  (Ar^_.^  is  obtained  from  Eq. 
(19)  in  terms  of  the  stored  APj  vectors,  i  =  1,...,n).  Furthermore, 
the  rate  of  convergence  is  governed  by  the  condition  number  of  matrix 
A  instead  of  A^A,  which  results  in  a  substantial  improvement  over  the 
normal  equations.  Again,  as  in  the  normal  equations,  the  rate  of 
convergence  and  number  of  iterations  is  significantly  improved  by 
preconditioning  Eq.  (10)  by  the  diagonal  elements  and  making  good 
extrapolations  for  from  previous  frequencies. 


The  rate  of  convergence  can  slow  down  in  the  asymmetric  algorithm 

when  the  angle  between  and  AP^  is  sufficiently  small  that  the 

incremental  reductions  in  the  residue  from  successive  iterations  is 

negligible.  If  this  condition  ever  occurs  before  desired  convergence  is 

reached,  the  iteration  loop  is  restarted  with  the  latest  approximation  as 

the  initial  estimate  and  the  initial  search  direction  in  Eq.  (19)  is 

modified  for  the  restart  to  maximize  (AP  ,  r  -).  One  method  to 

n  n-1 
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modify  the  search  directions  is  to  replace  r^  .j  in  Eqs.  (19)  and  (20) 
with  A^r^_^,  which  essentially  switches  the  method  to  the  symmetric 
case  for  one  iteration.  Again,  this  procedure  is  still  in  testing  and 
other  possibilities  may  be  equally  plausible. 


'.-V-J 


•  *  ^  ’J 
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?.3  COMPUTATIONAL  COMPARISONS 


The  remainder  of  this  section  will  discuss  why  BINTEQ  represents  the 

optimal  BIE  formulation  in  terms  of  execution  time  to  generate  the  exact 

solution.  A  theoretically  straightforward  inversion  to  solve  the  linear 

system  in  Eq.  (2)  would  have  been  extremely  computational  inefficient 

and  numerically  ill-conditioned  because  it  would  have  entailed  solving  an 

enormous  singular  system  of  (6N+3)q  x  (6N+3)q  complex  equations  for 

the  (6N+3)q  unknown  displacement  and  traction  boundary  values,  with 

"N"  being  the  number  of  layers  and  "q"  being  the  average  number  of 

nodes  on  one  interface.  Any  full  inversion  type  approach  would 

3  3 

require  on  the  order  of  (6N+3)  q  floating  point  operations  per 

2 

frequency  compared  to  about  144Nq  P  for  the  BINTEQ  solution 

technique,  with  P  being  the  number  of  iterations.  For  example,  with 

1 1 

N=5  and  q=256,  the  full  inversion  requires  about  5x10  operations  per 

9 

frequency  whereas  the  BINTEQ  technique  would  require  about  1x10 

operations  with  P=20  iterations  per  frequency.  Furthermore,  the  overall 

cost  of  a  BIfJTEQ  calculation  is  controlled  by  the  speed  of  repetitively 

multiplying  3qx3q  complex  matrices  times  3qx1  complex  vectors  and  is 

ideally  suited  to  execute  with  an  array  processor.  Using  an  FPS 

AP-120B  array  processor,  this  sample  problem  would  take  about  100 

seconds  per  frequency,  using  BINTEQ,  \.ihich  is  significantly  faster 

than  a  typical  elastic  finite  difference  type  calculation  with  full  volume 

gridding.  Compared  to  other  BIE  techniques  and  assuming  adaptability 

to  achieve  the  high  compute  rates  of  the  referenced  array  processor, 

Ferguson's  full  inversion  algorithm  (1982)  is  estimated  to  take  about  13 

hours  per  frequency  for  this  problem;  Schuster's  approach  (1984)  is 

estimated  to  take  about  4  hours  per  frequency  for  this  problem  based 
3  2 

on  (1728N+432)q  +  72Nq  P  floating  point  operations  per  frequency 

because  of  the  N  complex  6qx6q  interface  interaction  matrix  inversions 
and  the  complex  3qx3q  free  surface  interaction  matrix  inversion. 


3.0  RESULTS 


3.1  VALIDATION  CALCULATIONS 


Since  the  improved  method  is  still  in  the  final  development  and  testing 
stages,  the  validation  exercises  from  the  original  BINTEQ  method  will  be 
described  in  this  section.  There  are  very  few  exact  solutions  available 
for  checking  the  algorithm.  Therefore,  the  philosophy  of  the  validation 
stage  has  been  to  perform  as  many  internal  checks  as  possible  while  at 
the  same  time  comparing  solutions  against  those  from  established  and 
available  algorithms. 

An  exhaustive  set  of  internal  checks  of  the  BINTEQ  algorithm  have 
been  completed  for  the  2-D  acoustic  code:  verification  that  the  free 
surface  boundary  concitions  and  the  continuity  conditions  are  being 
satisfied  for  a  wide  range  of  problems;  verification  that  the  effects  of 
varying  impedance  contrasts  from  no  contrast  to  a  range  of  contrasts 
are  correct;  verification  that  the  lower  limit  of  two  nodal  points  per 
wavelength  produces  the  same  results  as  finer  sampling  for  the  same 
structure  (although  user  may  occasionally  use  finer  sampling  to  avoid 
spatial  aliasing  with  structures  having  very  steeply  sloping 
irregularities);  verification  that  there  are  no  spurious  edge  reflections 
contaminating  solution  for  a  wide  range  of  conditions  including  very 
high  material  attenuation  (Q  =  10)  to  almost  no  material  attenuation  (Q  = 
2000);  verification  that  extra  padding  at  edge  of  model  has  no  effect  on 
results;  and  numerous  other  verifications. 

External  validations  have  included:  verification  of  all  arrivals  for 

problems  with  exactly  Mat  interfaces  including  head  waves  and  Stoneley 
waves  against  Sierra  Geophysics’  VESPA-p|y|  algorithm  (Apsel,  1982) 
which  simulates  complete  solutions  in  flat  multilayered  viscoelastic 
structures;  verification  of  geometrical  arrival  times  and  amplitudes  for  a 
wide  range  of  problems  with  irregular  interfaces  against  Sierra 
Geophysics'  OUIK.^1^  three-dimensional  raytracing  algorithm  (Lundquist, 
el  al.,  198?);  and  vjrification  of  all  arrivals  including  diffractions 
against  a  limited  set  of  available  finite  difference  calculations. 


30 


Figure  3  shows  the  model  geometry  for  the  flat  layer  comparison  to  the 
exact  VESPA  solution.  The  source  is  located  just  beneath  the  interface 
at  a  depth  of  350  meters  at  an  offset  of  500  meters.  Receivers  are 
located  across  the  free  surface  from  offset*  of  600  out  to  2100  meters. 
The  compressional  velocities  are  5  km/sec  and  10  km/sec  for  layers  1 
and  2;  densities  are  2.0  gm/cc  and  2.5  gn/cc;  quality  factors  (Q)  are 
100  and  1000;  the  layer  thickness  is  300  me'ers.  The  simulations  include 
frequencies  from  10  to  100  Hz.  Figures  4  and  5  show  the  BINTEQ  and 
VESPA  shot  records,  respectively  for  all  80  receivers  and  travel  times 
ranging  from  0.0  to  0.8  seconds.  The  traces  have  been  convolved  with 
a  Ricker  wavelet  with  a  center  frequency  of  50  Hz.  The  BINTEQ 
solution  in  Figure  4  has  been  scaled  by  a  time  ramp  to  amplify  the  late 
arrivals  and  the  3-D  VESPA  solution  in  Figure  5  has  been  scaled  by  an 
extra  square  root  of  time  to  approximately  account  for  the  differences 
in  2-D  and  3-D  geometrical  spreading.  The  first  event  is  the  direct 
compressional  arrival  that  is  refracted  up  to  the  free  surface  through 
the  interface.  The  event  moving  out  at  the  slower  5  km/sec  velocity 
starting  at  0.1  seconds  and  group  21  is  a  non-geometrical  arrival  that 
"tunnels"  through  the  interface.  The  fact  that  BINTEQ  is  able  to  model 
such  non-geomelrical  arrivals  is  quite  encouraging  for  the  method.  The 
next  four  events  beginning  at  0.18,  0.30,  0.42  and  0.54  are  multiples  of 
order  1  through  4  reverberating  in  the  top  layer.  The  excellent 
agreement  with  the  exact  VESPA  solution  lends  considerable  confidence 
in  the  method. 

Figures  6,  7,  and  8  repeat  the  same  validation  study  except  that  a  200 
meter  thick  basin  has  been  placed  on  the  interface  as  shown  in  Figure 
6.  This  time  the  BINTEQ  solution  in  Figure  7  is  being  compared 
against  the  QUIK  raytracing  solution  in  Figure  8  which  provides  only 
geometrical  arrivals.  All  other  parameters  are  the  same  as  for  the  flat 
layer  comparison.  The  agreement  is  excellent  for  the  geometrical 
arrivals  modeled  by  the  QUIK  raytracing  code.  Notice  how  BINTEQ 
gels  continuous  events  including  all  detracted  energy  off  the  corners  of 
the  basin  whereas  the  raytracing  code  is  unable  to  model  the  diffracted 
energy.  For  example,  the  gap  in  the  QUIK  solution  at  0.1  seconds  and 
groups  29  through  36  is  where  the  smaller  diffracted  event  through  the 
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FLAT  LAYER  MODEL  FOR  BINTEQA^ESPA  COMPARISON 


Free 
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Figure  3 

Model  geometry  for  the  flat  layer  validation  study  between  BINTEQ  and  VESPA  shown  in  Figures  4  and  5. 
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GROUP  NUHBER 

Rgur*  4 

BINTEQ  soliiOon  for  tfw  flat  modei  shown  in  Figim  3. 


VESPR:  COMPRRISON  TO  BINTEQ  FOR  SIMPLE  2-LnYER  MODEL 

HORIZONTRL  SEISMIC  SECTION;  DEPTH  =  20.00 

SZZ  FOR  EXPLOSION  DEPTH  =  350.0 
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VESPA  Mact  solution  for  ths  flat  mocial  shown  in  Rgura  3. 


BASIN  MODEL  FOR  BINTEQ/QUIK  COMPARISON 


Figure  6 

Model  geometry  for  the  basin  model  comparison  between  BINTEQ  and  QUIK  shown  in  Figures  7  and  8. 
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GROUP  NUMBER 

Figure  7 

BINTEQ  solution  for  the  basin  model  shown  In  Figure  6. 


:SHOT:  SEISMIC  SHOT  RECORD  FOR  BRSIN  MODEE 

GEOMETRIC  nRRIVRLS  FOR  NEHR-SURFRCE  GROUPS 
PRESSURE  COMPONENT  FOR  POINT  SOURCE  RT  350  METER  DEPTH 


GROUP  NUMBER 

Figure  8 

QUIK  raytradng  solution  for  the  basin  model  shown  in  Figure  6. 


left  corner  of  the  basin  appears.  The  non-geometrical  event  that 
arrives  early  at  the  largest  offsets  in  the  BINTEQ  solution  is  a  head 
wave  traveling  at  the  faster  velocity  of  10  km/sec  along  the  bottom  of 
the  basin.  It  is  also  interesting  to  note  all  the  back  scattered  energy 
bouncing  around  in  the  basin. 

The  model  for  the  last  set  of  validation  studies  shown  in  this  section  is 
presented  at  the  top  of  Figure  9.  This  more  complex  low-velocity 
wedge  model  was  used  by  Koslov  and  Baysal  (1982)  to  demonstrate  the 
accuracy  of  the  finite  difference  technique  against  a  physical  model.  In 
the  physical  model,  the  low  velocity  wedge  was  submerged  in  water  and 
pinched  out  against  a  plexiglass  plate.  The  scaled  seismic  velocities 
used  in  the  comparison  are  4  km/sec  in  the  top  layer  (water),  2.3 
km/sec  in  the  second  layer  (wedge)  and  6  km/sec  in  the  underlying 
layer  (plate).  The  source  is  located  under  the  free  surface  at  a  depth 
of  100  meters  and  an  offset  of  4440  meters.  Receivers  are  located  just 
under  the  free  surface  at  a  depth  of  50  meters  at  offsets  ranging 
between  750  to  5700  meters  at  an  increment  of  100  meters.  The  BINTEQ 
shot  record  is  shown  at  the  same  scale  as  the  model  in  Figure  9  and 
includes  all  arrivals  in  the  frequency  range  of  4  to  30  Hz.  The 
synthetic  seismograms  are  shown  as  a  function  of  time  from  0.25  to  4.75 
seconds  so  as  to  clip  off  the  targe  impulse  for  the  receivers  near  the 
source  at  zero  time.  The  most  important  events  are  identified  by 
numbers  1  through  17. 

Event  1  is  simply  the  direct  arrival  from  the  source  plus  the  reflection 
off  the  free  surface.  Events  2  through  6  represent  arrivals  interacting 
with  the  wedge  and  the  free  surface:  events  2,  4  and  6  are  the 
primary  reflection  and  the  first  and  second  multiple  reflections  off  the 
wedge,  respectively;  event  3  represents  the  diffracted  arrivals  off 
corner  B;  and  event  5  represents  a  twice  diffracted  arrival  off  corner 
B. 

Events  7  through  12  represent  arrivals  involving  primary  reflections  off 
the  plate;  event  7  is  the  primary  off  the  plate  transmitting  through 
and  diffracting  off  the  wedge  in  the  vicinity  of  corner  C  and  traveling 


resolve  the  sharp  corner  of  the  wedge  with  too  few  samples.  The 
comparison  was  repeated  with  the  4th  order  finite  difference  algorithm 
of  John  Vidale  at  California  Institute  of  Technology.  The  results  from 
this  comparison  are  shown  in  Figure  13  on  the  same  scale  as  the 
BINTEQ  and  QUIK  shot  records  in  Figures  11  and  12.  The  agreement 
with  this  finite  difference  simulation  is  quite  spectacular. 
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BINTEO<  SCISHIC  SHOT  RECORD  FOR  LOM-VELOCITY  HEDGE  MODEL 
COHPLETe  ACOUSTIC  SOLUTION  FOR  NEAR-SURFflCE  CROUPS 
PRESSURE  COMPONENT  FOR  EXPLOSION  DEPTH  OF  100  METERS 


RCE 


FOURIER  FINITE-DIFFERENCE  METHOD 


FINITE-DIFFERENCE  METHOD 

SHOT  RECORD  FOR  LOW-VELOCITY  WEDGE  MODEL 


GROUP  NUMBER 

Figure  13 

Finite  difference  solution  of  Vidale  and  Clayton  (1984)  for  the  wedge  model  shown  in  Figure  9. 


3.2  CALCULATION  FOR  AFGL  MODELS 


The  results  from  two  sample  calculations  are  presented  in  Figures  14 
and  15  to  show  the  effects  of  surface  irregularities  on  the  synthetic 
seismograms.  Both  2-D  models  represent  simple  layer  over  half-space 
structures  with  a  flat  free  surface  and  an  irregular  interface  with 

compressional  velocities  of  4  and  8  km/sec,  densities  of  2.0  and  3.4 
gm/cc  and  material  quality  fa  Tors  of  50  and  1000  in  the  layer  and  the 
underlying  half-space,  respectively.  The  simulations  have  been  carried 
out  from  0  to  2.5  Hz  and  include  all  possible  arrivals  in  the  time 

window  of  0  to  50  seconds.  The  models  are  shown  to  the  left  of  the 
seismic  sections  with  tick  roarks  on  the  free  surface  denoting  the 
location  of  receivers  at  model  distances  from  8  km  to  94.4  km.  The 
interfaces  have  been  discretized  at  a  sampling  interval  of  0.8  km  for 
128  points  satisfying  the  minimum  requirement  of  2  sample  points  per 
wavelength  at  the  highest  frequency  of  interest.  Both  figures  show 
true  relative  amplitudes  except  for  a  multiplicative  scaling  factor  of  the 
hypocentral  distance  to  the  power  unity  for  display  purposes  only. 

The  results  in  Figure  14  can  be  used  to  examine  the  effects  of 

irregularities  on  refracting  horizons.  A  point  isotropic  source  is 
located  at  a  horizontal  model  distance  of  10  km  at  a  depth  of  11  km, 
just  one  kilometer  below  the  irregular  interface.  The  results  will  be 
discussed  with  the  aide  of  the  labelled  arrivals.  First  one  can  follow 
the  direct  arrival  (A)  as  it  moves  out  at  the  upper  velocity  of  4 
km/sec.  At  about  20  km,  .irrival  B  emerges  ahead  of  arrival  A  as  it 
travels  almost  horizontally  at  the  higher  velocity  of  8  km/sec  before 
being  refracted  up  to  the  free  surface.  \s  arrival  B  passes  through 
the  interface  irregularities,  it  undergoes  significant  interference  effects 
destroying  the  coherence  of  the  phase  and  causing  some  back  scattered 
energy  (arrivals  C).  Then  after  passing  through  the  irregularities,  a 
coherent  arrival  re-emerges  with  a  moveout  of  8  km/sec  (arrival  D). 
The  first  multiple  in  the  layer  is  arrival  E  and  it  is  interesting  to  watch 
it  merge  with  the  direct  arrival  A  at  distances  beyond  about  70  km. 
Arrivals  F  and  I  represent  back-scattered  energy  from  the  first  multiple 
off  the  interface  irregularities.  At  distances  beyond  about  72  km,  the 
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Figure  15.  BINTEQ  simulation  (0-2.5  hz)  of  the  seismic  section 
for  a  point  source  located  in  the  model  shown  at  the  left  for 
a  deep  basin  structure. 
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first  multiple  (arrival  E)  generates  a  head  wave  on  the  interface 
(arrival  G)  with  a  moveout  of  8  km/sec.  This  arrival  appears  to  grow 
beyond  90  km  because  of  constructive  interference  from  direct  arrivals 
scattered  toward  the  free  surface  at  a  moveout  of  4  km/sec  (arrivals 
originating  between  C  and  D).  Arrivals  H  and  K  represent  the  next  two 
higher  order  multiples  in  the  layer  and  arrivals  J  plus  other  unmarked 
small  arrivals  represent  more  back  scattered  energy  off  the  interface 
irregularities  from  these  higher  multiples. 

The  results  in  Figure  15  can  be  used  to  understand  the  effects  of  wave 
propagation  through  a  basin.  Direct  arrivals  A  and  B  have  the  same 
interpretation  as  in  Figure  14.  In  this  figure,  arrival  C  represents  a 
true  head  wave  off  the  bottom  of  the  basin  and  arrival  D  represents  a 
creeping  wave  corresponding  to  the  continuation  of  arrival  D  past  the 
eastern  rise  of  the  basin.  The  sequence  of  arrivals  E  represent 
diffractions  off  the  rough  corners  of  the  discretely  sampled  western 
slope  of  the  basin  (and  are  not  present  at  lower  frequency  where 
sampling  is  sufficient  to  emulate  a  smooth  slope).  As  in  the  previous 
example,  arrival  F  represents  back-scattered  energy  off  the  eastern 
rise  of  the  basin.  Arrivals  G,  H  and  I  are  the  first  multiple,  its 
back-scattered  energy  off  eastern  slope  and  its  head  wave  generated 
after  the  eastern  slope,  analogous  to  arrivals  E,  F  and  G  of  the 
previous  example.  Arrivals  J  and  K  represent  the  next  two  higher 
order  multiples  and  again  the  small  back-scattered  arrivals  from  these 
multiples  are  too  small  to  demarcate. 

In  conclusion,  BINTEQ  has  been  verified  through  an  exhaustive  series 
of  internal  and  external  validation  tests.  Exact  synthetic  seismograms 
were  shown  for  two  models  to  study  the  influence  of  interface 
irregularities  on  seismic  wave  propagation  in  the  earth.  The  most 
important  findings  from  these  studies  are;  1)  the  interference  effects 
observed  from  interface  irregularities  (arrivals  C  in  Figure  2)  may  help 
explain  the  phase  incoherence  observed  in  real  data  from  head  wave 
energy  along  the  Moho  discontinuity;  (2)  significant  back-scattered 
energy  is  to  be  expected  from  interface  irregularities  (arrivals  F,  I,  J 
in  Figure  14)  and  from  basins  (arrivals  F,  H  in  Figure  15);  and  (3)  the 
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appearance  and  disappearance  of  critically  refracted  energy  (arrivals  G 
in  Figure  14  and  I  in  Figure  15)  could  easily  have  been  misinterpreted 
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